Let us consider the dataframe wages, which containes information about 3294 USA working individuals. The data are taken from the National Longitudinal Survey and are related to 1987. The variable as are listed below and the output of the str command is given
## 'data.frame': 3296 obs. of 5 variables:
## $ exper : int 9 12 11 9 8 9 8 10 12 7 ...
## $ male : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...
## $ school: int 13 12 11 14 14 14 12 12 10 12 ...
## $ wage : num 6.32 5.48 3.64 4.59 2.42 ...
## $ region: Factor w/ 3 levels "Center","North",..: 1 1 3 1 1 1 1 1 1 3 ...
The aim of the study is to analyze the potential relationship between the response variable wage and the explanatory variables considered in the dataframe. Describe how to perform a preliminary data analysis on this dataframe, using suitable R commands. Moreover, consider the following plots and discuss the possibility of measuring the variable wage in the logarithmic scale
source("./../functions.R", local = knitr::knit_global())
summary(wages)
## exper male school wage
## Min. : 1.000 female:1569 Min. : 3.00 Min. : 0.07656
## 1st Qu.: 7.000 male :1727 1st Qu.:11.00 1st Qu.: 3.62673
## Median : 8.000 Median :12.00 Median : 5.20578
## Mean : 8.042 Mean :11.63 Mean : 5.81639
## 3rd Qu.: 9.000 3rd Qu.:12.00 3rd Qu.: 7.30723
## Max. :18.000 Max. :16.00 Max. :112.79193
## region
## Center:1164
## North :1105
## South :1027
##
##
##
numeric.graph(wages, c(4,1,3), qqplot = 4, no.density = c(1,3), draw.normal.curve = c(1,3))
## [1] "Skewness indexes"
## wage exper school
## 9.5491015 0.1555681 -0.3934935
## [1] "Kurtosis indexes - 3"
## wage exper school
## 212.9197746 0.6528604 1.2394911
categorical.graph(wages, c(2,5))
La variabile risposta presenta una elevatissima simmetria e curtosi, con i valori praticamente tutti concentrati a sinistra del grafico con la presenza di due outlier elevatissimi. Il qqplot conferma la presenza di una coda destra molto pesante, mentre i quantili centrali seguono abbastanza bene la distribuzione normale. Si consiglia perciò o di rimuovere gli outlier più lontani oppure di introdurre una trasformata.
exper invece sembra seguire piuttosto bene la curva della distribuzione normale (buona simmetria, curtosi leggermente ipernormale).
Considerazioni simili valgono anche per school (anche qui buona simmetria), la distribuzione di presenta con una frequenza del valore centrale (11) molto elevata e una coda sinistra leggermente appesantita.
Vista la natura di queste due variabili non ha molto senso considerare gli outlier.
wages$logwage <- log(wages$wage)
numeric.graph(wages, 6, qqplot = 6)
## [1] "Skewness indexes"
## [1] -1.11844
## [1] "Kurtosis indexes - 3"
## [1] 4.28369
La trasformata logaritmica migliora notevolmente sia simmetria che curtosi della variabile, tuttavia appesantisce molto la coda sinistra (infatti la campana rimane fortemente ipernormale), per cui eventuali assunzioni di normalità sulla variabile vanno fatte con cautela.
boxplots.graph(6, c(1:3,5), wages)
Dai grafici non sembrano esserci delle particolari correlazioni tra
log(wage) e le altre variabili.
boxplots.graph(1, 2, wages)
mean(wages$exper[wages$male == 'male'])
## [1] 8.323104
mean(wages$exper[wages$male == 'female'])
## [1] 7.732314
Dal grafico e dal confronto delle medie sembra esserci una (debole) correlazione tra exper e male.
#Rimuovo outlier
newWage <- wages[wages$wage < 30,]
newWage$logwage <- log(newWage$wage)
newWage$wageT <- transformResponseWithBoxCox(newWage$wage~1, newWage, 4)
## [1] "Ideal lambda for wage : 0.4"
numeric.graph(newWage, c(6,7), qqplot = c(6,7))
## [1] "Skewness indexes"
## logwage wageT
## -1.23982911 0.04971715
## [1] "Kurtosis indexes - 3"
## logwage wageT
## 4.1493448 0.8880249
La trasformata di Box Cox con \(\labda=0.4\) sembra preferibile rispetto a quella logaritmica.
In order to describe the effect of the factor male on the response log(wage) we may analyze this plot, where the probability distribution of log(wage) is represented by considering the kernel density estimates conditioned on the two levels (1 male, 0 female) of the variable male
Il grafico conferma l’assenza di correlazione tra log(wage) e male, infatti le due curve di densità sono praticamente sovrapponibili.
With the commands mod.0<-lm(log(wage) ∼ male,data=wages) and mod.1<-lm(log(wage) ∼ exper*male, data=wages), two regression models are defined for describing the potential effect of male and exper on the response log(wage). Comment the model fitting outcomes given by the function summary (Hint: consider the fact that the average years of experience in the sample is lower for women than for men).
mod.0 <- lm(log(wage) ~ male,data=wages)
summary(mod.0)
##
## Call:
## lm(formula = log(wage) ~ male, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0445 -0.3073 0.0544 0.3839 3.0325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.47475 0.01559 94.59 <2e-16 ***
## malemale 0.21826 0.02154 10.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6176 on 3294 degrees of freedom
## Multiple R-squared: 0.03023, Adjusted R-squared: 0.02994
## F-statistic: 102.7 on 1 and 3294 DF, p-value: < 2.2e-16
mod.1 <- lm(log(wage) ~ exper*male, data=wages)
summary(mod.1)
##
## Call:
## lm(formula = log(wage) ~ exper * male, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0906 -0.3050 0.0560 0.3792 3.0468
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.193551 0.057470 20.768 < 2e-16 ***
## exper 0.036367 0.007156 5.082 3.94e-07 ***
## malemale 0.463707 0.079062 5.865 4.93e-09 ***
## exper:malemale -0.032071 0.009518 -3.369 0.000762 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6153 on 3292 degrees of freedom
## Multiple R-squared: 0.03792, Adjusted R-squared: 0.03704
## F-statistic: 43.25 on 3 and 3292 DF, p-value: < 2.2e-16
Il confronto tra i due modelli mostra come l’interazione tra exper e male (rappresentato dalla variabile dummy malemale) sia rilevante nel modello, di conseguenza lo sono anche i loro main effects.
Finally, a complete regression model is fitted mod.2 <- lm(log(wage) ∼., data=wages) and the following output is obtained by the R function summary.
summary(lm(log(wage) ~. -logwage, data=wages))
##
## Call:
## lm(formula = log(wage) ~ . - logwage, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0008 -0.2821 0.0468 0.3673 3.2337
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.279709 0.090321 -3.097 0.00197 **
## exper 0.034618 0.004549 7.610 3.55e-14 ***
## malemale 0.246474 0.020607 11.961 < 2e-16 ***
## school 0.122909 0.006278 19.578 < 2e-16 ***
## regionNorth 0.051107 0.024505 2.086 0.03709 *
## regionSouth 0.047168 0.024969 1.889 0.05898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5831 on 3290 degrees of freedom
## Multiple R-squared: 0.1364, Adjusted R-squared: 0.1351
## F-statistic: 103.9 on 5 and 3290 DF, p-value: < 2.2e-16
Describe how to interpret these results, and then suggest how to proceed with further analyses.
Il modello qui presentato mostra come tutti le variabili siano regressori rilevanti per logwage.
plot(lm(log(wage) ~. -logwage, data = wages))
Il modello presenta alcune criticità dal punto di vista della normalità dei residui. Le ipotesi di omoschedasticità e linearità sono rispettate.
Cosa otteniamo utilizzando sempre la trasformata logaritmica senza outliers
fitLog <- lm(logwage~exper+male+school+region, data = newWage)
summary(fitLog)
##
## Call:
## lm(formula = logwage ~ exper + male + school + region, data = newWage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9969 -0.2792 0.0490 0.3710 1.8645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.299749 0.089186 -3.361 0.000786 ***
## exper 0.036006 0.004493 8.013 1.54e-15 ***
## malemale 0.243299 0.020352 11.954 < 2e-16 ***
## school 0.123402 0.006196 19.915 < 2e-16 ***
## regionNorth 0.056040 0.024193 2.316 0.020597 *
## regionSouth 0.045451 0.024669 1.842 0.065500 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5754 on 3285 degrees of freedom
## Multiple R-squared: 0.1405, Adjusted R-squared: 0.1392
## F-statistic: 107.4 on 5 and 3285 DF, p-value: < 2.2e-16
plot(fitLog)
drop1(lm(wageT~exper+male+school+region, data = newWage), test = "F")
## Single term deletions
##
## Model:
## wageT ~ exper + male + school + region
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 3300.7 21.66
## exper 1 55.62 3356.3 74.66 55.3587 1.275e-13 ***
## male 1 170.70 3471.4 185.60 169.8906 < 2.2e-16 ***
## school 1 444.86 3745.5 435.76 442.7430 < 2.2e-16 ***
## region 2 3.29 3304.0 20.94 1.6363 0.1949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dal risultato di drop1 risulta che per la nuova trasformata il regressore region non è significativo.
fitBoxCox <- lm(wageT~exper+male+school, data=newWage)
summary(fitBoxCox)
##
## Call:
## lm(formula = wageT ~ exper + male + school, data = newWage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7995 -0.6243 -0.0142 0.6259 4.4220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.997868 0.153070 -6.519 8.16e-11 ***
## exper 0.057779 0.007824 7.384 1.93e-13 ***
## malemale 0.462388 0.035457 13.041 < 2e-16 ***
## school 0.226935 0.010796 21.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.003 on 3287 degrees of freedom
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.1516
## F-statistic: 196.9 on 3 and 3287 DF, p-value: < 2.2e-16
plot(fitBoxCox)
La trasformata calcolata con il metodo di Box-Cox permette di rientrare (seppur parzialmente) nell’ipotesi di normalità dei residui.
AIC(fitLog, fitBoxCox)
## df AIC
## fitLog 7 5710.026
## fitBoxCox 5 9362.390
BIC(fitLog, fitBoxCox)
## df BIC
## fitLog 7 5752.719
## fitBoxCox 5 9392.884
Gli indici \(R^2\) sono più alti nella trasformata Box Cox rispetto al logaritmo, tuttavia le statistiche AIC e BIC tendono a preferire quest’ultima.